Cite this Article
Zhao Lu-Si, Chen Chun-Ping, Liu Lin-Lin, Yu Hong-Xia, Chen Yi, Wang Xiao-Chun. Magnetism and piezoelectricity of hexagonal boron nitride with triangular vacancy. Chinese Physics B, 2018, 27(1): 016301  
Permissions
Magnetism and piezoelectricity of hexagonal boron nitride with triangular vacancy
Zhao Lu-Si1, 2, Chen Chun-Ping1, 2, Liu Lin-Lin1, 2, Yu Hong-Xia1, 2, Chen Yi1, 2, Wang Xiao-Chun1, 2, †
Institute of Atomic and Molecular Physics, Jilin University, Changchun 130012, China
Jilin Provincial Key Laboratory of Applied Atomic and Molecular Spectroscopy (Jilin University), Changchun 130012, China

 

† Corresponding author. E-mail: wangxiaochun@jlu.edu.cn

Abstract

First-principle calculations reveal that the configuration system of hexagonal boron nitride (h-BN) monolayer with triangular vacancy can induce obvious magnetism, contrary to that of the nonmagnetic pristine boron nitride monolayer. Interestingly, the h-BN with boron atom vacancy (VB-BN) displays metallic behavior with a total magnetic moment being 0.46μB per cell, while the h-BN with nitrogen atom vacancy (VN-BN) presents a half-metallic characteristic with a total magnetic moment being 1.0μB per cell. Remarkably, piezoelectric stress coefficient e11 of the VN-BN is about 1.5 times larger than that of pristine h-BN. Furthermore, piezoelectric strain coefficient d11 (12.42 pm/V) of the VN-BN is 20 times larger than that of pristine h-BN and also one order of magnitude larger than the value for the h-MoS2 monolayer, which is mainly due to the spin-down electronic state in the VN-BN system. Our study demonstrates that the nitrogen atom vacancies can be an efficient route to tailoring the magnetic and piezoelectric properties of h-BN monolayer, which have promising performances for potential applications in nano-electromechanical systems (NEMS) and nanoscale electronics devices.

PACS: 63.20.dk;75.70.Ak;77.65.Ly;85.85.+j
Keyword:first-principle;magnetism;piezoelectric;nano-electromechanical systems
1. Introduction

The successful synthesis of atomically thick graphene.[1] has triggered enormous interest in research on low-dimensional nanomaterials due to its exciting unusual properties and promising applications in material science. Consequently, hexagonal boron nitride (h-BN) monolayer, a two-dimensional material with the same honeycomb structure as graphene, exhibits completely different physical properties from graphene[2–4] and has been brought to the forefront of research.[5–11] The advances in experimental synthesis of h-BN sheet have led us to explore the properties of h-BN sheet by decorating its surface. High-resolution transmission electron microscopy studies have revealed that the defects in h-BN sheets are primarily triangle-shaped vacancy defects caused during the fabrication of h-BN.[12]

Actually, the structural and electronic properties,[13,14] chemical reactivity[15] and molecular adsorption and dissociation behaviors[16,17] of vacancy-decorated h-BN have been reported. In the pristine h-BN sheet, the charge transfer from B to N atoms and the orbital hybridization make electrons paired and thus the system is nonmagnetic.[18] Interestingly, the vacancies in the h-BN sheet can induce interesting magnetic properties.[19–21] However, how the magnetic property of h-BN sheet changes by increasing the density of the vacancies and considering the interaction between the neighboring vacancies is still an unsolved problem.

The rapidly growing demand for nano-electromechanical systems (NEMS) and nanoscale electronics devices calls for nanoscale piezoelectric materials, promoting investigations into low-dimensional materials.[22–24] For a crystal to be piezoelectric, it must have no center of inversion symmetry.[25] Graphene lacks any intrinsic piezoelectricity due to its symmetry inversion center. None of many two-dimensional (2D) monolayer structures exhibits an inversion symmetry, and thus all of them are intrinsically piezoelectric.[26–29] The most prominent 2D structures are hexagonal boron nitride (h-BN) and some transition metal dichalcogenides (TMDCS) due to their large piezoelectric coefficients.[26,28,29] Interestingly, it was found that that the modified h-BN or h-AlN sheet by hydrogen and/or fluorine also shows an enhanced piezoelectric performance by first-principles calculation.[30,31] The ab initio calculations indicate that the emergence of piezoelectricity in g-C3N4 is due to the fact that a stable phase of graphene nitride nanosheet is riddled with regularly spaced triangular holes.[32] We predict that the h-BN sheet with triangle-shaped holes may improve piezoelectric performance. Electronic and magnetic properties of h-BN monolayer with mono-vacancy defects have been reported frequently.[13,20,21,33] However, so far, magnetic and piezoelectric properties of h-BN riddled with triangle-shaped holes have not been systematically investigated yet.

Here, we use first-principle calculation to investigate the magnetic and piezoelectric properties of the vacancy-decorated h-BN monolayer with an array of densest vacancies which have noncentrosymmetric triangular holes in structure. These decorated structures include VB-BN and VN-BN, which indicate the h-BN monolayer with one boron atom and one nitrogen atom vacancy decorated in one super cell, respectively. These h-BN monolayer systems with different triangle-shaped vacancies could potentially introduce new functionalities through cross-coupling effects between electrical, mechanical and magnetic responses.

2. Computational details

The first-principle calculations are based on spin-polarized density functional theory using projector-augmented wave methods[34] as implemented in the Vienna ab initio simulation package (VASP).[35–37] Exchange–correlation interactions are described in the generalized gradient approximation (GGA)[38] within the Perdew–Burke–Ernzerhof (PBE).[39] Plane waves with a kinetic energy cutoff of 600 eV are used to expand the valance electron (2s22p1 for B, 2s22p3 for N) wavefunctions. A vacuum spacing of 15 Å is used to avoid the interactions between the layers. The Brillouin zone is sampled by a 7 × 7 × 1 k-point grid within the Monkhorst–Pack scheme[40] for our calculations, which ensures that the total energy convergence is within 0.01% difference with respect to the number of k points. For geometry relaxation, the total energy change and the force on each atom are less than 10−5 eV and 10−3 eV/Å, respectively. The elastic tensor coefficients (Cijkl) are calculated by using the finite difference method[41] and the strain coefficients of the piezoelectric tensor (eijk) are calculated by using density-functional perturbation theory (DFPT).[41]

3. Results and discussion
3.1. Structural and stable properties

Following the idea that the h-BN monolayer is decorated by the triangle-shaped vacancy, we discuss the two configuration systems as shown in Fig. 1. Compared with the pristine h-BN sheet, the h-BN with boron vacancies (marked as VB-BN) consists of those vacancies each of which is due to the missing of one boron atom in primitive cell as shown in Fig. 1(a). As a result, in view of the overall sheet, this forms a regular spaced triangle-shaped vacancy. The three doubly coordinated nitrogen atoms are the closest ones to the vacancy. For convenience of discussion, we here denote the trebly coordinated and the doubly coordinated N atoms as N1 and N2, respectively. The unit cell of the VB-BN contains seven atoms as denoted by black dashed lines in Fig. 1(a) in which the B-to-N ratio is 3:4. The optimized lattice constant is a = 4.975Å. The B–N1 (1.403Å) and B–N2 (1.408Å) bond lengths of the VB-BN are less than the B–N (1.45Å) bond length of the pristine h-BN sheet. This fact indicates that the interactions between the B and N atoms of the VB-BN are stronger than that of the pristine h-BN. Here, we perform a Bader charge analysis for better elucidating the effective charges at the B and N atoms. The different outer shell electronic configurations and electronegativities of atoms cause an imbalance of electric charges within the B–N bond. Consequently, the electrons are shared unevenly, forming a polar covalent bond. A positive number denoted by “+” means that this atom loses electrons and has a positive charge, and similarly, a negative number denoted by “–” means that this atom obtains electrons and has a negative charge. As listed in Fig. 1(a), Bader analysis reveals that there are certain electronic charge transfers from the less electronegative B atoms to the more electronegative N atoms, which is the same as the pristine h-BN. In the unit cell of the VB-BN, about 2.18|e| charge is transferred to the N1 and N2 atoms from each B atom. Each N1 atom obtain 2.15|e| while the N2 atom gains less charge (1.46 ± 0.01|e|).

Fig. 1. (color online) Optimized structures of (a) boron vacancy VB-BN and (b) nitrogen vacancy VN-BN. Gray and green spheres represent nitrogen and boron atoms, respectively. Rhombus marked by black dashed lines denotes primitive cell.

Figure 1(b) displays the structure model of the h-BN with nitrogen vacancies (denoted as VN-BN). Like the VB-BN defect case, the VN-BN defect loses one nitrogen atom in the unit cell compared with the pristine h-BN sheet as shown in Fig. 1(b). In contrast to the VB-BN, it is clear that the nitrogen vacancy is surrounded by the three doubly coordinated B atoms. Geometrical relaxation of the VN-BN is performed by considering the rhombic primitive cell as denoted by black dashed lines in Fig. 1(b). Here, we group the trebly coordinated and the doubly coordinated B atoms as B1 and B2, respectively. The optimized lattice constant of the VN-BN is found to be a = 4.940Å, which is smaller than that of the VB-BN. In the unit cell of the VN-BN, the ratio between the numbers of atoms B and N is 4:3. The B1–N bond length is 1.419Å which is slightly larger than the B–N1 (1.403Å) bond length of the VB-BN, while B2-N bond length is 1.456Å which is close to the B–N (1.450Å) bond length of the pristine h-BN sheet. Therefore, it seems that the VB-BN is easier to form than the VN-BN. The Bader charge analysis demonstrates that each N atom gains 2.18|e| from the neighboring B atoms. The B1 atom transfers 2.16|e| to the neighboring N atoms, but the B2 atom transfers 1.47 or 1.46|e| to the neighboring N atoms.

3.2. Electronic and magnetic properties

Figure 2(a) shows the spin-polarized band structure and total density of states (TDOS) for the VB-BN system. We find that the bands of the spin-down electron states exceed slightly those of the spin-up electron states below an energy of 2 eV in the spin-polarized band structure, inducing a spontaneous polarization which is in accordance with the TDOS of this system. Unlike the semiconducting property of pristine h-BN sheet, this VB-BN system displays metallic behavior, which originates from the pseudodangling bonds of N atoms neighboring with B atom vacancy. Whether its state is a spin-up state or spin-down state, the valence band crosses the Fermi level between all high symmetry points, spreading the whole Brillouin zone. According to the partial density of states (PDOS) in Fig. 2(b), we further shows that the magnetic property is mainly contributed by unsaturated 2p orbital electrons of N2 atoms, which is consistent with the results of 3D spin charge density distribution in Fig. 2(d).

Fig. 2. (color online) ((a), (b)) Spin-polarized band structures with a Fermi level sets at zero energy, total density of states (TDOS) and partial density of states (PDOS) of the VB-BN. (c) 2D plot of charge density difference of VB-BN. (d) 3D plot of spin charge density distribution of the VB-BN atan isovalue of 0.004 e/Å3.

Three-dimensional (3D) plots of spin charge density distributions of the VB-BN in Fig. 2(d) visualize the magnetic property clearly. The 3D spin charge density distribution [Δρ(r) = ρ↑(r) − ρ↓(r)] reveals that the magnetic moments are mainly localized on N2 atoms and secondarily localized on N1. The total magnetic moment is 0.46μB per cell. The N1 and N2 atoms contribute an atomic moment of 0.059μB, 0.100μB, respectively. Compared with N1 and N2 atoms, each B atom carry very small magnetic moment (−0.011μB). In order to understand the chemical bonding of the VB-BN system, we draw its 2D plots of the charge density difference (defined as the total charge density minus the contributions of isolated B and N atoms) in Fig. 2(c). A positive value (solid contour line) indicates an increase in electron density, and a negative value (dashed contour line) means the electron density loss with respect to the superposition of the atomic electron density. Figure 2(c) shows that there is electronic charge accumulation in the middle of the bond between the neighboring B and N atoms, which indicates how the B–N covalent bonds are formed. There is electronic charge depletion around the N2 atoms. Because each nitrogen atom has a valence configuration of 2s22p3, two electrons of the N2 atom participate in the covalent bond with the boron atoms, two electrons of the N2 atom make up a lone pair electron, and the remaining one electron is shared between pseudodangling sp3 bond which leads each N2 atom to own a magnetic moment. Because three electrons of the N1 atom participate in the covalent bond with the boron atoms, there is no electron to be shared between pseudodangling sp3 bonds for the N1 atom. However, the shared electron of N2 atom has an influence on the homologous N1 atom and thus causes the N1 atom to have a small magnetic moment.

The spin-polarized band structure and total density of states (TDOS) of the VN-BN are plotted in Fig. 3(a). In the spin-polarized band structure, the spin-up and spin-down states have different dispersions. The bands with energy lower than −4 eV are fully occupied and thus contribute little to spin polarization. However, the flat band near Fermi energy is split into two branches. The spin-up branch is occupied and the spin-down branch is empty, leading to a spontaneous polarization with a total magnetic moment of 1.0μB per cell, for which the B1 and B2 atoms contribute atomic moments of 0.072μB and 0.132μB, respectively. Comparing with the VB-BN, the magnetic moment of the VN-BN is about twice larger than that of the VB-BN, but the leading role in their magnetism is played by these doubly coordinated atoms surrounding the vacancy for both VB-BN and VN-BN systems. Interestingly, the VN-BN presents a half-metallic characteristic where spin-down state remains as an indirect band-gap insulator with a band gap of about 4.2 eV, which is smaller than 4.68 eV for h-BN by 0.48 eV.[26] The spin-up state is metallic as shown in the spin-polarized band structure. The result is in accordance with the calculated TDOS. In the spin-down state, the TDOS is zero at the Fermi level with an energy gap of 4.2 eV, while there is a finite TDOS at the Fermi level in the spin-up state, revealing half-metallic behavior in the VN-BN system. Thus, the VN-BN is only conductive to the spin-up electrons, which leads to a 100% spin polarization. The corresponding spin charge density distribution with an isosurface of 0.01 e/Å3 is depicted in Fig. 3(d), which shows that the spin states at different atoms have different spin orientations and the magnetism is mainly contributed by the B atoms. The 2D plots of the charge density difference of the VN-BN system in Fig. 3(c) depicts the nature of the chemical bonding. Figure 3(c) shows that the B atom surrounded by three triangular N vacancies has stronger interaction with neighboring N atoms than other B atoms. This fact can contribute to the stability of the VN-BN structure. As is well known, B atom has a valence configuration of 2s22p1, once one nitrogen atom is removed, each B2 atom will remain one unpaired electron. Therefore the unpaired electrons form pseudodangling bond at the position of vacancy and the B2 atoms will produce exactly identical magnetic moment. Like N1 atom in the VB-BN, the B1 atom in VN-BN has a little magnetic moment which is affected by B2 atoms. The point is demonstrated by the calculated partial density of states in Fig. 3(b) again, which are mainly dominated by 2p states of B2 atom at the Fermi level.

Fig. 3. (color online) The legend is similar to that in Fig. 2 while the results in panels (a)–(d) are for the VN-BN. (d) Yellow and light blue indicate the spin-up and spin-down states, respectively, in the isosurface 0.01 e/Å3.
3.3. Mechanical and piezoelectric properties

In preceding discuss, we find the VB-BN to be metallic and hence not piezoelectric. Then we perform the piezoelectric effect calculation on VN-BN system. Adopting these optimized structures, we can calculate the elastic constants and piezoelectric coefficients of the pristine and N atom vacancy-decorated h-BN. The elastic constants reflect the stress–strain relation of material, which is defined as

where σ and ε are the stress and strain tensor, respectively, and Cij are the elastic constants.

Similarly, the piezoelectric coefficients can be expressed as

where both the piezoelectric stress coefficients eijk and the piezoelectric strain coefficients dijk are third-rank tensors, relating the polarizations Pi to strain εjk and stress σjk, respectively.

We know that the performing of symmetry analysis on the elastic and piezoelectric tensors can further reduce the number of independent tensor coefficients by using the point group of the 2D material. There are many regular triangle-shaped holes in vacancy-decorated h-BN. However, compared with the pristine h-BN, the vacancy-decorated h-BN still belongs to the D3h () point group. Employing the standard Voigt notation, i.e., 1-xx, 2-yy, 3-zz, 4-yz, 5-xz, 6-xy, the elastic and piezoelectric coefficients can become

where we can conclude that the only unique piezoelectric coefficients are e11 and d11. In addition, the corresponding piezoelectric strain coefficients d11 are further obtained from the inversion of the transformation relationship as follows:

The elastic constants and piezoelectric coefficients calculated for the pristine h-BN and the VN-BN are listed in Table 1. In the present study, the elastic constants and piezoelectric coefficients of the pristine h-BN sheet are calculated to be C11 = 289 N/m and C12 = 64 N/m by using the finite differences method and e11 = 1.40 × 10−10 C/m and d11 = 0.62 pm/V by using density-functional perturbation theory (DFPT), which are well consistent with previous results of C11 = 291 N/m, C12 = 62 N/m, e11 = 1.3810−10 C/m, and d11 = 0.60 pm/V by Duerloo et al.[26] Obviously, the elastic constants of the VN-BN become smaller than those of the pristine counterpart. Comparing the C11 of the VN-BN with that of the pristine h-BN, we find that the elastic constant 289 N/m of the pristine h-BN is much larger than triple the elastic constant (86 N/m) of the VN-BN, indicating that the induced N atom vacancies make it more flexible. As can be seen from Table 1, the calculated piezoelectric stress coefficients e11 of the VN-BN monolayer are of the same order of magnitude as other popular piezoelectric materials listed for comparison. The VN-BN exhibits a greater piezoelectric response than the pristine h-BN It is seen that the piezoelectric stress coefficient e11 of VN-BN (2.36 × 10−10 C/m) is much larger than that of the pristine BN monolayer (1.40 × 10−10 C/m). It is worth noting that the piezoelectric stress coefficient e11 of VN-BN is larger than 1.5 times that of pristine h-BN, larger than the maximum one predicted for F- and H-doped h-BN in a chair conformation.[31] Furthermore, the piezoelectric stress coefficient e11 of VN-BN is larger than those of 2H-NbSe2 (2.22 × 10−10 C/m) and 2H-NbTe2 (1.84 × 10−10 C/m). The calculated piezoelectric strain coefficient d11 of VN-BN (12.42 pm/V) is one order of magnitude higher than those of all the other materials listed in Table 1, and it is 20 times larger than that of pristine h-BN and also 1 order of magnitude higher than the values for the h-MoS2, 2H-NbSe2 and 2H-NbTe2 monolayer. Figure 3(a) shows that the spin-up state of the VN-BN presents metallic characteristic, and the spin-down state of the VN-BN remains an indirect band-gap insulator with a band gap of about 4.2 eV. The piezoelectric coefficient is mainly due to the electronic structure of spin-down electronic states in the VN-BN system. The regularly triangular vacancies do not change the symmetry of the h-BN sheet when the symmetric center is located at the atom in the lower-left corner of the unit cell. Then the VN-BN still belongs to the D3h () point group. However, the non-centrosymmetric triangular vacancy will strengthen the polarized degree of this sheet. That is why the piezoelectric stress coefficients e11 of the VN-BN is larger than that of the pristine h-BN sheet. We can deduce that other structures like the h-BN with losing one atom to form triangular vacancy in the unit cell could make larger piezoelectric response. Moreover, the pristine h-BN and MoS2 monolayer lose their piezoelectric effect when even-numbered layers are present, but the piezoelectricities of these h-BN monolayers riddled with regularly spaced triangular holes, which are studied here and listed in Table 1, can be retained if they are stacked like g-C3N4 based on a similar mechanism.[32]

Table 1.

Calculated elastic constants (Cij) and piezoelectric coefficients (eij and dij) of the pristine and VN-BN. Other experimental and theoretical values of popular 2D piezoelectric materials are listed for comparison. Elastic constants (Cij) are in units of N/m, piezoelectric stress coefficients (eij) in 10−10 C/m for 2D system, and piezoelectric strain coefficients (dij) in pm/V.

.
4. Conclusions

In this work, we investigate the magnetic and piezoelectric properties of hexagonal boron nitride (h-BN) monolayer riddled with regularly spaced triangular holes by using first-principle calculations. It is found that the triangular holes can indeed induce spontaneous magnetizations in the VB-BN and VN-BN systems. More interestingly, the band structure of VN-BN indicates that the VN-BN has the desired half-metallic behavior. As for piezoelectricity, the piezoelectric strain coefficients d11 (12.42 pm/V) of the VN-BN is also 20 times larger than that of pristine h-BN and one order of magnitude higher than the value of the h-MoS2 monolayer. The spin-down electronic states mainly contribute to the higher piezoelectric coefficient in the VN-BN system. It is worth noting that the VN-BN system riddled with different regularly spaced triangular holes has no center of inversion symmetry and may not lose piezoelectricity when even-numbered layers are present. Our study demonstrates that the triangle-shaped vacancy is used to tailor the magnetic and piezoelectric properties of h-BN monolayer. The co-existence of piezoelectricity and magnetism in h-BN with triangle-shaped vacancy means promising performance for potential applications in the integration of magnetic and electromechanical nano-devices.

Reference
[1] Novoselov K S Geim A K Morozov S V Jiang D Zhang Y Dubonos S V Grigorieva I V Firsov A A 2004 Science 306 666
[2] Watanabe K Taniguchi T Kanda H 2004 Nat. Mater. 3 404
[3] Novoselov K Jiang D Schedin F Booth T Khotkevich V Morozov S Geim A 2005 Proc. Natl. Acad. Sci. USA 102 10451
[4] Degheidy A R Elkenany E B 2017 Chin. Phys. 26 086103
[5] Han W Q Wu L Zhu Y Watanabe K Taniguchi T 2008 Appl. Phys. Lett. 93 223103
[6] Nag A Raidongia K Hembram K P Datta R Waghmare U V Rao C 2010 ACS Nano 4 1539
[7] Zhi C Bando Y Tang C Kuwahara H Golberg D 2009 Adv. Mater. 21 2889
[8] Jasuja K 2011 Designing Nanoscale Constructs from Atomic Thin Sheets of Graphene, Boron Nitride and Gold Nanoparticles for Advanced Material Applications Kansas State University
[9] Zhao L S Wang Y Chen C P Liu L L Yu H X Zhang Y Chen Y Wang X C 2017 Physica 91 82
[10] Liu L L Chen C P Zhao L S Wang Y Wang X C 2017 Carbon 115 773
[11] Chen X Xu L Liu L L Zhao L S Chen C P Zhang Y Wang X C 2017 Appl. Surf. Sci. 396 1020
[12] Jin C Lin F Suenaga K Iijima S 2009 Phys. Rev. Lett. 102 195505
[13] Azevedo S Kaschny J R de Castilho C M C de Brito Mota F 2009 The European Physical Journal 67 507
[14] Azevedo S Kaschny J R de Castilho C M de Brito Mota F 2007 Nanotechnology 18 495707
[15] Lin Y Williams T V Cao W Elsayed-Ali H E Connell J W 2010 The Journal of Physical Chemistry 114 17434
[16] Noorizadeh S Shakerzadeh E 2012 Computational Materials Science 56 122
[17] Lvova N Ananina O Y 2016 Comput. Mater. Sci. 115 11
[18] Zhou J Wang Q Sun Q Jena P 2010 Phys. Rev. 81 085442
[19] Yin J Li J Hang Y Yu J Tai G Li X Zhang Z Guo W 2016 Small 12 2942
[20] Wang Y Ding Y 2014 J. Phys.: Condens. Matt. 26 43
[21] Si M S Xue D S 2007 Phys. Rev. 75 193409
[22] Ekinci K Roukes M 2005 Rev. Sci. Instrum 76 061101
[23] Huang X M H Zorman C A Mehregany M Roukes M L 2003 Nature 421 496
[24] Sai N Mele E 2003 Phys. Rev. 68 241405
[25] Reed E J 2015 Nat. Nanotechnology 10 106
[26] Duerloo K A N Ong M T Reed E J 2012 The Journal of Physical Chemistry Letters 3 2871
[27] Michel K Verberck B 2009 Phys. Rev. 80 224301
[28] Li W Li J 2015 Nano Research 8 3796
[29] Wang J L Hu W D 2017 Chin. Phys. 26 037106
[30] Ding Y Wang Y 2016 J. Mater. Chem. 4 1517
[31] Noor-A-Alam M Kim H J Shin Y H 2014 Phys. Chem. Chem. Phys. 16 6575
[32] Zelisko M Hanlumyuang Y Yang S Liu Y Lei C Li J Ajayan P M Sharma P 2014 Nat. Commun. 5 4284
[33] Yang J Kim D Hong J Qian X 2010 Surf. Sci. 604 1603
[34] Blöchl P E 1994 Phys. Rev. 50 17953
[35] Kresse G Hafner J 1994 Phys. Rev. 49 14251
[36] Kresse G Furthmüller J 1996 Phys. Rev. 54 11169
[37] Kresse G Furthmüller J 1996 Comput. Mater. Sci. 6 15
[38] Perdew J P Chevary J A Vosko S H Jackson K A Pederson M R Singh D J Fiolhais C 1992 Phys. Rev. 46 6671
[39] Perdew J P Burke K Ernzerhof M 1996 Phys. Rev. Lett. 77 3865
[40] Monkhorst H J Pack J D 1976 Phys. Rev. 13 5188
[41] Wu X Vanderbilt D Hamann D 2005 Phys. Rev. 72 035105
[42] Blonsky M N Zhuang H L Singh A K Hennig R G 2015 ACS Nano 9 9885
[43] Zhu H Y Wang Y Xiao J Liu M Xiong S M Wong Z J Ye Z L Ye Y Yin X B Zhang X 2015 Nat. Nanotechnol. 10 151